function [like,grad]=consumpsearchstrucmlogit(b,Y,lambda,ut,grad_4yr,gpdiff,AdjNG,AdjG,sdemog,S,q,Beta)

    flag = Y>0;

    in_white_collar = ismember(Y,[2 4 7 9 12 14 17 19]);
    prev_WC         = ut.sprevs(:,end);

    b2flg   = 2+[1:ut.number2];
    b4sflg  = 2+[1+ut.number2:ut.number2+ut.number4s];
    b4nsflg = 2+[1+ut.number2+ut.number4s:ut.number2+ut.number4s+ut.number4ns];
    bwptflg = 2+[1+ut.number2+ut.number4s+ut.number4ns:ut.number2+ut.number4s+ut.number4ns+ut.numberpt];
    bwftflg = 2+[1+ut.number2+ut.number4s+ut.number4ns+ut.numberpt:ut.number2+ut.number4s+ut.number4ns+ut.numberpt+ut.numberft];
    bwcflg  = 2+[1+ut.number2+ut.number4s+ut.number4ns+ut.numberpt+ut.numberft:ut.number2+ut.number4s+ut.number4ns+ut.numberpt+ut.numberft+ut.numberwc];

    alpha     = b(1);
    rho       = b(2);
    b2        = b(b2flg);
    b2temp    = [b2(1:sdemog+1);0;b2(sdemog+2:end-7);0;0;b2(end-6:end)];% two 0s: grad_4yr, white_collar
    b2tflg    = find(b2temp~=0);
    b4s       = b(b4sflg);
    b4stemp   = [b4s(1:sdemog+1);0;b4s(sdemog+2:end-7);0;0;b4s(end-6:end)];
    b4stflg   = find(b4stemp~=0);
    b4ns      = b(b4nsflg);
    b4nstemp  = [b4ns(1:sdemog+1);0;b4ns(sdemog+2:end-7);0;0;b4ns(end-6:end)];
    b4nstflg  = find(b4nstemp~=0);
    bwpt      = b(bwptflg);
    bwpttemp  = [bwpt(1:sdemog);0;0;bwpt(sdemog+1:end-3);0;0;0;0;bwpt(end-2:end)];% only four 0s here bc need to estimate white_collar parameter
    bwpttflg  = find(bwpttemp~=0);
    bwft      = b(bwftflg);
    bwfttemp  = [bwft(1:sdemog);0;0;bwft(sdemog+1:end-3);0;0;0;0;0;bwft(end-2:end)];
    bwfttflg  = find(bwfttemp~=0);
    bwc       = b(bwcflg);
    bwctemp   = [bwc(1:sdemog);0;0;bwc(sdemog+1:end-3);0;0;0;0;0;bwc(end-2:end)];
    bwctflg   = find(bwctemp~=0);

    % indices to flag consumption and non-consumption covariates in matrices
    cidx           = sdemog+2;

    % "utility" for 2-year college:
    vng2ftbc = (ut.X2ftbc(flag,:)*(b2temp+bwfttemp        )+ut.X2ftbc(flag,cidx)*alpha)+((1-Beta).*gpdiff(flag,1)*rho)+AdjNG(flag,1);
    vng2ftwc = (ut.X2ftwc(flag,:)*(b2temp+bwfttemp+bwctemp)+ut.X2ftwc(flag,cidx)*alpha)+((1-Beta).*gpdiff(flag,2)*rho)+AdjNG(flag,2);
    vng2ptbc = (ut.X2ptbc(flag,:)*(b2temp+bwpttemp        )+ut.X2ptbc(flag,cidx)*alpha)+((1-Beta).*gpdiff(flag,3)*rho)+AdjNG(flag,3);
    vng2ptwc = (ut.X2ptwc(flag,:)*(b2temp+bwpttemp+bwctemp)+ut.X2ptwc(flag,cidx)*alpha)+((1-Beta).*gpdiff(flag,4)*rho)+AdjNG(flag,4);
    vng2     = (ut.X2nw(  flag,:)*(b2temp                 )+ut.X2nw(  flag,cidx)*alpha)+((1-Beta).*gpdiff(flag,5)*rho)+AdjNG(flag,5);

    % "utility" for 4-year college: Science majors
    vng4sftbc = (ut.X4sftbc(flag,:)*(b4stemp+bwfttemp        )+ut.X4sftbc(flag,cidx)*alpha)+((1-Beta).*gpdiff(flag,6 )*rho)+AdjNG(flag,6);
    vng4sftwc = (ut.X4sftwc(flag,:)*(b4stemp+bwfttemp+bwctemp)+ut.X4sftwc(flag,cidx)*alpha)+((1-Beta).*gpdiff(flag,7 )*rho)+AdjNG(flag,7);
    vng4sptbc = (ut.X4sptbc(flag,:)*(b4stemp+bwpttemp        )+ut.X4sptbc(flag,cidx)*alpha)+((1-Beta).*gpdiff(flag,8 )*rho)+AdjNG(flag,8);
    vng4sptwc = (ut.X4sptwc(flag,:)*(b4stemp+bwpttemp+bwctemp)+ut.X4sptwc(flag,cidx)*alpha)+((1-Beta).*gpdiff(flag,9 )*rho)+AdjNG(flag,9);
    vng4s     = (ut.X4snw(  flag,:)*(b4stemp                 )+ut.X4snw(  flag,cidx)*alpha)+((1-Beta).*gpdiff(flag,10)*rho)+AdjNG(flag,10);

    % Non-Science majors
    vng4nsftbc = (ut.X4nsftbc(flag,:)*(b4nstemp+bwfttemp        )+ut.X4nsftbc(flag,cidx)*alpha)+((1-Beta).*gpdiff(flag,11)*rho)+AdjNG(flag,11);
    vng4nsftwc = (ut.X4nsftwc(flag,:)*(b4nstemp+bwfttemp+bwctemp)+ut.X4nsftwc(flag,cidx)*alpha)+((1-Beta).*gpdiff(flag,12)*rho)+AdjNG(flag,12);
    vng4nsptbc = (ut.X4nsptbc(flag,:)*(b4nstemp+bwpttemp        )+ut.X4nsptbc(flag,cidx)*alpha)+((1-Beta).*gpdiff(flag,13)*rho)+AdjNG(flag,13);
    vng4nsptwc = (ut.X4nsptwc(flag,:)*(b4nstemp+bwpttemp+bwctemp)+ut.X4nsptwc(flag,cidx)*alpha)+((1-Beta).*gpdiff(flag,14)*rho)+AdjNG(flag,14);
    vng4ns     = (ut.X4nsnw(  flag,:)*(b4nstemp                 )+ut.X4nsnw(  flag,cidx)*alpha)+((1-Beta).*gpdiff(flag,15)*rho)+AdjNG(flag,15);

    % "utility" for working
    vngwptbc = (ut.Xngwptbc(flag,:)*(bwpttemp        )+ut.Xngwptbc(flag,cidx)*alpha)+((1-Beta).*gpdiff(flag,16)*rho)+AdjNG(flag,16);
    vngwptwc = (ut.Xngwptwc(flag,:)*(bwpttemp+bwctemp)+ut.Xngwptwc(flag,cidx)*alpha)+((1-Beta).*gpdiff(flag,17)*rho)+AdjNG(flag,17);
    vngwftbc = (ut.Xngwftbc(flag,:)*(bwfttemp        )+ut.Xngwftbc(flag,cidx)*alpha)+((1-Beta).*gpdiff(flag,18)*rho)+AdjNG(flag,18);
    vngwftwc = (ut.Xngwftwc(flag,:)*(bwfttemp+bwctemp)+ut.Xngwftwc(flag,cidx)*alpha)+((1-Beta).*gpdiff(flag,19)*rho)+AdjNG(flag,19);

    % "utility" for grad school 
    vgwptbc = (ut.Xgwptbc( flag,:)*(bwpttemp                )+ut.Xgwptbc( flag,cidx)*alpha)+(0*rho)+AdjG(flag,16);
    vgwptwc = (ut.Xgwptwc( flag,:)*(bwpttemp+bwctemp        )+ut.Xgwptwc( flag,cidx)*alpha)+(0*rho)+AdjG(flag,17);
    vgwftbc = (ut.Xgwftbc( flag,:)*(bwfttemp                )+ut.Xgwftbc( flag,cidx)*alpha)+(0*rho)+AdjG(flag,18);
    vgwftwc = (ut.Xgwftwc( flag,:)*(bwfttemp+bwctemp        )+ut.Xgwftwc( flag,cidx)*alpha)+(0*rho)+AdjG(flag,19);

    %%% debugging:
    %allmatNG = [vng2ftbc vng2ftwc vng2ptbc vng2ptwc vng2 vng4sftbc vng4sftwc vng4sptbc vng4sptwc vng4s vng4nsftbc vng4nsftwc vng4nsptbc vng4nsptwc vng4ns vngwptbc vngwptwc vngwftbc vngwftwc];
    %allmatG = [vgwptbc vgwptwc vgwftbc vgwftwc];
    %disp('sum stats X4nsftbc');
    %sumopt = struct('Weights',q(flag));
    %summarize(ut.X4nsftbc(flag,:),sumopt);
    %disp('sum stats X4nsftwc');
    %sumopt = struct('Weights',q(flag));
    %summarize(ut.X4nsftwc(flag,:),sumopt);
    %disp('sum stats X4nsptbc');
    %sumopt = struct('Weights',q(flag));
    %summarize(ut.X4nsptbc(flag,:),sumopt);
    %disp('sum stats X4nsptwc');
    %sumopt = struct('Weights',q(flag));
    %summarize(ut.X4nsptwc(flag,:),sumopt);
    %disp('sum stats X4ns');
    %sumopt = struct('Weights',q(flag));
    %summarize(ut.X4nsnw(flag,:),sumopt);
    %disp('sum stats non-grad valfn');
    %sumopt = struct('Weights',q(flag).*(grad_4yr(flag)==0));
    %summarize(allmatNG,sumopt);
    %disp('sum stats grad valfn');
    %sumopt = struct('Weights',q(flag).*(grad_4yr(flag)==1));
    %summarize(allmatG,sumopt);
    %disp('sum stats all valfn');
    %sumopt = struct('Weights',q(flag));
    %summarize([allmatNG allmatG],sumopt);
    %disp('sum stats prevWC');
    %sumopt = struct('Weights',q(flag));
    %summarize(prev_WC(flag),sumopt);
    %disp('sum stats non-grad Adj');
    %sumopt = struct('Weights',q(flag).*(grad_4yr(flag)==0));
    %summarize(AdjNG(flag,1:19),sumopt);
    %disp('sum stats grad Adj');
    %sumopt = struct('Weights',q(flag).*(grad_4yr(flag)==1));
    %summarize(AdjG(flag,16:19),sumopt);
    %disp('choice parameters');
    %b
    %blablabla

    %% Form likelihood (3 parts)
    % Part 1a
    log_dem_ngno = log(1+exp(vng2ftbc)+exp(vng2ptbc)+exp(vng2)+exp(vng4sftbc)+exp(vng4sptbc)+exp(vng4s)+exp(vng4nsftbc)+exp(vng4nsptbc)+exp(vng4ns)+exp(vngwptbc)+exp(vngwftbc));
    log_p2ftbc_ngno   = vng2ftbc  -log_dem_ngno;
    log_p2ptbc_ngno   = vng2ptbc  -log_dem_ngno;
    log_p2_ngno       = vng2      -log_dem_ngno;
    log_p4sftbc_ngno  = vng4sftbc -log_dem_ngno;
    log_p4sptbc_ngno  = vng4sptbc -log_dem_ngno;
    log_p4s_ngno      = vng4s     -log_dem_ngno;
    log_p4nsftbc_ngno = vng4nsftbc-log_dem_ngno;
    log_p4nsptbc_ngno = vng4nsptbc-log_dem_ngno;
    log_p4ns_ngno     = vng4ns    -log_dem_ngno;
    log_pwptbc_ngno   = vngwptbc  -log_dem_ngno;
    log_pwftbc_ngno   = vngwftbc  -log_dem_ngno;
    log_ph_ngno       =           -log_dem_ngno;

    % Part 1b
    log_dem_ngof = log(1+exp(vng2ftbc)+exp(vng2ftwc)+exp(vng2ptbc)+exp(vng2ptwc)+exp(vng2)+exp(vng4sftbc)+exp(vng4sftwc)+exp(vng4sptbc)+exp(vng4sptwc)+exp(vng4s)+exp(vng4nsftbc)+exp(vng4nsftwc)+exp(vng4nsptbc)+exp(vng4nsptwc)+exp(vng4ns)+exp(vngwptbc)+exp(vngwptwc)+exp(vngwftbc)+exp(vngwftwc));
    log_p2ftbc_ngof   = vng2ftbc  -log_dem_ngof;
    log_p2ftwc_ngof   = vng2ftwc  -log_dem_ngof;
    log_p2ptbc_ngof   = vng2ptbc  -log_dem_ngof;
    log_p2ptwc_ngof   = vng2ptwc  -log_dem_ngof;
    log_p2_ngof       = vng2      -log_dem_ngof;
    log_p4sftbc_ngof  = vng4sftbc -log_dem_ngof;
    log_p4sftwc_ngof  = vng4sftwc -log_dem_ngof;
    log_p4sptbc_ngof  = vng4sptbc -log_dem_ngof;
    log_p4sptwc_ngof  = vng4sptwc -log_dem_ngof;
    log_p4s_ngof      = vng4s     -log_dem_ngof;
    log_p4nsftbc_ngof = vng4nsftbc-log_dem_ngof;
    log_p4nsftwc_ngof = vng4nsftwc-log_dem_ngof;
    log_p4nsptbc_ngof = vng4nsptbc-log_dem_ngof;
    log_p4nsptwc_ngof = vng4nsptwc-log_dem_ngof;
    log_p4ns_ngof     = vng4ns    -log_dem_ngof;
    log_pwptbc_ngof   = vngwptbc  -log_dem_ngof;
    log_pwptwc_ngof   = vngwptwc  -log_dem_ngof;
    log_pwftbc_ngof   = vngwftbc  -log_dem_ngof;
    log_pwftwc_ngof   = vngwftwc  -log_dem_ngof;
    log_ph_ngof       =           -log_dem_ngof;

    % Part 2a
    log_dem_grno     = log(1+exp(vgwptbc)+exp(vgwftbc));
    log_pwptbc_grno  = vgwptbc-log_dem_grno;
    log_pwftbc_grno  = vgwftbc-log_dem_grno;
    log_ph_grno      =        -log_dem_grno;

    % Part 2b
    log_dem_grof     = log(1+exp(vgwptbc)+exp(vgwptwc)+exp(vgwftbc)+exp(vgwftwc));
    log_pwptbc_grof  = vgwptbc-log_dem_grof;
    log_pwptwc_grof  = vgwptwc-log_dem_grof;
    log_pwftbc_grof  = vgwftbc-log_dem_grof;
    log_pwftwc_grof  = vgwftwc-log_dem_grof;
    log_ph_grof      =        -log_dem_grof;

    log_like_ngno = (Y(flag)==1).*log_p2ftbc_ngno    +...
                    (Y(flag)==3).*log_p2ptbc_ngno    +...
                    (Y(flag)==5).*log_p2_ngno        +...
                    (Y(flag)==6).*log_p4sftbc_ngno   +...
                    (Y(flag)==8).*log_p4sptbc_ngno   +...
                    (Y(flag)==10).*log_p4s_ngno      +...
                    (Y(flag)==11).*log_p4nsftbc_ngno +...
                    (Y(flag)==13).*log_p4nsptbc_ngno +...
                    (Y(flag)==15).*log_p4ns_ngno     +...
                    (Y(flag)==16).*log_pwptbc_ngno   +...
                    (Y(flag)==18).*log_pwftbc_ngno   +...
                    (Y(flag)==20).*log_ph_ngno;

    log_like_ngof = (Y(flag)==1).*log_p2ftbc_ngof    +...
                    (Y(flag)==2).*log_p2ftwc_ngof    +...
                    (Y(flag)==3).*log_p2ptbc_ngof    +...
                    (Y(flag)==4).*log_p2ptwc_ngof    +...
                    (Y(flag)==5).*log_p2_ngof        +...
                    (Y(flag)==6).*log_p4sftbc_ngof   +...
                    (Y(flag)==7).*log_p4sftwc_ngof   +...
                    (Y(flag)==8).*log_p4sptbc_ngof   +...
                    (Y(flag)==9).*log_p4sptwc_ngof   +...
                    (Y(flag)==10).*log_p4s_ngof      +...
                    (Y(flag)==11).*log_p4nsftbc_ngof +...
                    (Y(flag)==12).*log_p4nsftwc_ngof +...
                    (Y(flag)==13).*log_p4nsptbc_ngof +...
                    (Y(flag)==14).*log_p4nsptwc_ngof +...
                    (Y(flag)==15).*log_p4ns_ngof     +...
                    (Y(flag)==16).*log_pwptbc_ngof   +...
                    (Y(flag)==17).*log_pwptwc_ngof   +...
                    (Y(flag)==18).*log_pwftbc_ngof   +...
                    (Y(flag)==19).*log_pwftwc_ngof   +...
                    (Y(flag)==20).*log_ph_ngof;

    log_like_grno = (Y(flag)==16).*log_pwptbc_grno  +...
                    (Y(flag)==18).*log_pwftbc_grno  +...
                    (Y(flag)==20).*log_ph_grno;

    log_like_grof = (Y(flag)==16).*log_pwptbc_grof  +...
                    (Y(flag)==17).*log_pwptwc_grof  +...
                    (Y(flag)==18).*log_pwftbc_grof  +...
                    (Y(flag)==19).*log_pwftwc_grof  +...
                    (Y(flag)==20).*log_ph_grof;

    like_o  = exp(log_like_ngof.*(grad_4yr(flag)==0)+log_like_grof.*(grad_4yr(flag)==1)); % equation G.4 in appendix
    like_no = exp(log_like_ngno.*(grad_4yr(flag)==0)+log_like_grno.*(grad_4yr(flag)==1)); % equation G.4 in appendix

    ll0     = log(lambda(flag).*like_o + (1-lambda(flag)).*like_no);
    ll1     = log(lambda(flag).*like_o);

    like = -ctranspose(q(flag))*( ( ll0.*(in_white_collar(flag)==0 & prev_WC(flag)==0) ) + ( ll1.*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) );


    %% Analytical Gradient
    dem_no  = lambda(flag).*like_o + (1-lambda(flag)).*like_no;

    % create matrix of terms that hold the q*(d-P)*grad terms that always show up in the gradient
    temp = zeros(sum(flag),25,4);
    temp(:,1 ,1) = (q(flag).*((Y(flag)==1 )-exp(log_p2ftbc_ngno  )).*(grad_4yr(flag)==0));
    temp(:,3 ,1) = (q(flag).*((Y(flag)==3 )-exp(log_p2ptbc_ngno  )).*(grad_4yr(flag)==0));
    temp(:,5 ,1) = (q(flag).*((Y(flag)==5 )-exp(log_p2_ngno      )).*(grad_4yr(flag)==0));
    temp(:,6 ,1) = (q(flag).*((Y(flag)==6 )-exp(log_p4sftbc_ngno )).*(grad_4yr(flag)==0));
    temp(:,8 ,1) = (q(flag).*((Y(flag)==8 )-exp(log_p4sptbc_ngno )).*(grad_4yr(flag)==0));
    temp(:,10,1) = (q(flag).*((Y(flag)==10)-exp(log_p4s_ngno     )).*(grad_4yr(flag)==0));
    temp(:,11,1) = (q(flag).*((Y(flag)==11)-exp(log_p4nsftbc_ngno)).*(grad_4yr(flag)==0));
    temp(:,13,1) = (q(flag).*((Y(flag)==13)-exp(log_p4nsptbc_ngno)).*(grad_4yr(flag)==0));
    temp(:,15,1) = (q(flag).*((Y(flag)==15)-exp(log_p4ns_ngno    )).*(grad_4yr(flag)==0));
    temp(:,16,1) = (q(flag).*((Y(flag)==16)-exp(log_pwptbc_ngno  )).*(grad_4yr(flag)==0));
    temp(:,18,1) = (q(flag).*((Y(flag)==18)-exp(log_pwftbc_ngno  )).*(grad_4yr(flag)==0));
    temp(:,20,1) = (q(flag).*((Y(flag)==20)-exp(log_ph_ngno      )).*(grad_4yr(flag)==0));

    temp(:,1 ,2) = (q(flag).*((Y(flag)==1 )-exp(log_p2ftbc_ngof  )).*(grad_4yr(flag)==0));
    temp(:,2 ,2) = (q(flag).*((Y(flag)==2 )-exp(log_p2ftwc_ngof  )).*(grad_4yr(flag)==0));
    temp(:,3 ,2) = (q(flag).*((Y(flag)==3 )-exp(log_p2ptbc_ngof  )).*(grad_4yr(flag)==0));
    temp(:,4 ,2) = (q(flag).*((Y(flag)==4 )-exp(log_p2ptwc_ngof  )).*(grad_4yr(flag)==0));
    temp(:,5 ,2) = (q(flag).*((Y(flag)==5 )-exp(log_p2_ngof      )).*(grad_4yr(flag)==0));
    temp(:,6 ,2) = (q(flag).*((Y(flag)==6 )-exp(log_p4sftbc_ngof )).*(grad_4yr(flag)==0));
    temp(:,7 ,2) = (q(flag).*((Y(flag)==7 )-exp(log_p4sftwc_ngof )).*(grad_4yr(flag)==0));
    temp(:,8 ,2) = (q(flag).*((Y(flag)==8 )-exp(log_p4sptbc_ngof )).*(grad_4yr(flag)==0));
    temp(:,9 ,2) = (q(flag).*((Y(flag)==9 )-exp(log_p4sptwc_ngof )).*(grad_4yr(flag)==0));
    temp(:,10,2) = (q(flag).*((Y(flag)==10)-exp(log_p4s_ngof     )).*(grad_4yr(flag)==0));
    temp(:,11,2) = (q(flag).*((Y(flag)==11)-exp(log_p4nsftbc_ngof)).*(grad_4yr(flag)==0));
    temp(:,12,2) = (q(flag).*((Y(flag)==12)-exp(log_p4nsftwc_ngof)).*(grad_4yr(flag)==0));
    temp(:,13,2) = (q(flag).*((Y(flag)==13)-exp(log_p4nsptbc_ngof)).*(grad_4yr(flag)==0));
    temp(:,14,2) = (q(flag).*((Y(flag)==14)-exp(log_p4nsptwc_ngof)).*(grad_4yr(flag)==0));
    temp(:,15,2) = (q(flag).*((Y(flag)==15)-exp(log_p4ns_ngof    )).*(grad_4yr(flag)==0));
    temp(:,16,2) = (q(flag).*((Y(flag)==16)-exp(log_pwptbc_ngof  )).*(grad_4yr(flag)==0));
    temp(:,17,2) = (q(flag).*((Y(flag)==17)-exp(log_pwptwc_ngof  )).*(grad_4yr(flag)==0));
    temp(:,18,2) = (q(flag).*((Y(flag)==18)-exp(log_pwftbc_ngof  )).*(grad_4yr(flag)==0));
    temp(:,19,2) = (q(flag).*((Y(flag)==19)-exp(log_pwftwc_ngof  )).*(grad_4yr(flag)==0));
    temp(:,20,2) = (q(flag).*((Y(flag)==20)-exp(log_ph_ngof      )).*(grad_4yr(flag)==0));

    temp(:,16,3) = (q(flag).*((Y(flag)==16)-exp(log_pwptbc_grno )).*(grad_4yr(flag)==1));
    temp(:,18,3) = (q(flag).*((Y(flag)==18)-exp(log_pwftbc_grno )).*(grad_4yr(flag)==1));
    temp(:,20,3) = (q(flag).*((Y(flag)==20)-exp(log_ph_grno     )).*(grad_4yr(flag)==1));

    temp(:,16,4) = (q(flag).*((Y(flag)==16)-exp(log_pwptbc_grof )).*(grad_4yr(flag)==1));
    temp(:,17,4) = (q(flag).*((Y(flag)==17)-exp(log_pwptwc_grof )).*(grad_4yr(flag)==1));
    temp(:,18,4) = (q(flag).*((Y(flag)==18)-exp(log_pwftbc_grof )).*(grad_4yr(flag)==1));
    temp(:,19,4) = (q(flag).*((Y(flag)==19)-exp(log_pwftwc_grof )).*(grad_4yr(flag)==1));
    temp(:,20,4) = (q(flag).*((Y(flag)==20)-exp(log_ph_grof     )).*(grad_4yr(flag)==1));

    grad_ngno=zeros(size(b));
    grad_ngno(1,1)       = -(ut.X2ftbc  (flag,cidx    ))'*( temp(:, 1,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no ) - ...
                            (ut.X2ptbc  (flag,cidx    ))'*( temp(:, 3,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X2nw    (flag,cidx    ))'*( temp(:, 5,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4sftbc (flag,cidx    ))'*( temp(:, 6,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4sptbc (flag,cidx    ))'*( temp(:, 8,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4snw   (flag,cidx    ))'*( temp(:,10,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4nsftbc(flag,cidx    ))'*( temp(:,11,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4nsptbc(flag,cidx    ))'*( temp(:,13,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4nsnw  (flag,cidx    ))'*( temp(:,15,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.Xngwptbc(flag,cidx    ))'*( temp(:,16,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.Xngwftbc(flag,cidx    ))'*( temp(:,18,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no);
    grad_ngno(2,1)       = -(     (1-Beta).*gpdiff(flag, 1      ))'*( temp(:, 1,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no ) - ...
                            (     (1-Beta).*gpdiff(flag, 3      ))'*( temp(:, 3,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (     (1-Beta).*gpdiff(flag, 5      ))'*( temp(:, 5,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (     (1-Beta).*gpdiff(flag, 6      ))'*( temp(:, 6,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (     (1-Beta).*gpdiff(flag, 8      ))'*( temp(:, 8,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (     (1-Beta).*gpdiff(flag,10      ))'*( temp(:,10,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (     (1-Beta).*gpdiff(flag,11      ))'*( temp(:,11,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (     (1-Beta).*gpdiff(flag,13      ))'*( temp(:,13,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (     (1-Beta).*gpdiff(flag,15      ))'*( temp(:,15,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (     (1-Beta).*gpdiff(flag,16      ))'*( temp(:,16,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (     (1-Beta).*gpdiff(flag,18      ))'*( temp(:,18,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no);
    grad_ngno(b2flg,1)   = -(ut.X2ftbc  (flag,b2tflg  ))'*( temp(:, 1,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X2ptbc  (flag,b2tflg  ))'*( temp(:, 3,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X2nw    (flag,b2tflg  ))'*( temp(:, 5,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no);
    grad_ngno(b4sflg,1)  = -(ut.X4sftbc (flag,b4stflg ))'*( temp(:, 6,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4sptbc (flag,b4stflg ))'*( temp(:, 8,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4snw   (flag,b4stflg ))'*( temp(:,10,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no);
    grad_ngno(b4nsflg,1) = -(ut.X4nsftbc(flag,b4nstflg))'*( temp(:,11,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4nsptbc(flag,b4nstflg))'*( temp(:,13,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4nsnw  (flag,b4nstflg))'*( temp(:,15,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no);
    grad_ngno(bwptflg,1) = -(ut.X2ptbc  (flag,bwpttflg))'*( temp(:, 3,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4sptbc (flag,bwpttflg))'*( temp(:, 8,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4nsptbc(flag,bwpttflg))'*( temp(:,13,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.Xngwptbc(flag,bwpttflg))'*( temp(:,16,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no);
    grad_ngno(bwftflg,1) = -(ut.X2ftbc  (flag,bwfttflg))'*( temp(:, 1,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4sftbc (flag,bwfttflg))'*( temp(:, 6,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4nsftbc(flag,bwfttflg))'*( temp(:,11,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.Xngwftbc(flag,bwfttflg))'*( temp(:,18,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no);

    grad_ngof_o=zeros(size(b));
    grad_ngof_o(1,1)       = -(ut.X2ftbc  (flag,cidx    ))'*( temp(:, 1,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X2ftwc  (flag,cidx    ))'*( temp(:, 2,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X2ptbc  (flag,cidx    ))'*( temp(:, 3,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X2ptwc  (flag,cidx    ))'*( temp(:, 4,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X2nw    (flag,cidx    ))'*( temp(:, 5,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4sftbc (flag,cidx    ))'*( temp(:, 6,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4sftwc (flag,cidx    ))'*( temp(:, 7,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4sptbc (flag,cidx    ))'*( temp(:, 8,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4sptwc (flag,cidx    ))'*( temp(:, 9,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4snw   (flag,cidx    ))'*( temp(:,10,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsftbc(flag,cidx    ))'*( temp(:,11,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsftwc(flag,cidx    ))'*( temp(:,12,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsptbc(flag,cidx    ))'*( temp(:,13,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsptwc(flag,cidx    ))'*( temp(:,14,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsnw  (flag,cidx    ))'*( temp(:,15,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.Xngwptbc(flag,cidx    ))'*( temp(:,16,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.Xngwptwc(flag,cidx    ))'*( temp(:,17,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.Xngwftbc(flag,cidx    ))'*( temp(:,18,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.Xngwftwc(flag,cidx    ))'*( temp(:,19,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) );
    grad_ngof_o(2,1)       = -(     (1-Beta).*gpdiff(flag, 1      ))'*( temp(:, 1,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (     (1-Beta).*gpdiff(flag, 2      ))'*( temp(:, 2,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (     (1-Beta).*gpdiff(flag, 3      ))'*( temp(:, 3,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (     (1-Beta).*gpdiff(flag, 4      ))'*( temp(:, 4,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (     (1-Beta).*gpdiff(flag, 5      ))'*( temp(:, 5,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (     (1-Beta).*gpdiff(flag, 6      ))'*( temp(:, 6,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (     (1-Beta).*gpdiff(flag, 7      ))'*( temp(:, 7,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (     (1-Beta).*gpdiff(flag, 8      ))'*( temp(:, 8,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (     (1-Beta).*gpdiff(flag, 9      ))'*( temp(:, 9,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (     (1-Beta).*gpdiff(flag,10      ))'*( temp(:,10,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (     (1-Beta).*gpdiff(flag,11      ))'*( temp(:,11,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (     (1-Beta).*gpdiff(flag,12      ))'*( temp(:,12,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (     (1-Beta).*gpdiff(flag,13      ))'*( temp(:,13,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (     (1-Beta).*gpdiff(flag,14      ))'*( temp(:,14,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (     (1-Beta).*gpdiff(flag,15      ))'*( temp(:,15,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (     (1-Beta).*gpdiff(flag,16      ))'*( temp(:,16,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (     (1-Beta).*gpdiff(flag,17      ))'*( temp(:,17,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (     (1-Beta).*gpdiff(flag,18      ))'*( temp(:,18,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (     (1-Beta).*gpdiff(flag,19      ))'*( temp(:,19,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) );
    grad_ngof_o(b2flg,1)   = -(ut.X2ftbc  (flag,b2tflg  ))'*( temp(:, 1,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X2ftwc  (flag,b2tflg  ))'*( temp(:, 2,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X2ptbc  (flag,b2tflg  ))'*( temp(:, 3,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X2ptwc  (flag,b2tflg  ))'*( temp(:, 4,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X2nw    (flag,b2tflg  ))'*( temp(:, 5,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) );
    grad_ngof_o(b4sflg,1)  = -(ut.X4sftbc (flag,b4stflg ))'*( temp(:, 6,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4sftwc (flag,b4stflg ))'*( temp(:, 7,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4sptbc (flag,b4stflg ))'*( temp(:, 8,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4sptwc (flag,b4stflg ))'*( temp(:, 9,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4snw   (flag,b4stflg ))'*( temp(:,10,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) );
    grad_ngof_o(b4nsflg,1) = -(ut.X4nsftbc(flag,b4nstflg))'*( temp(:,11,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsftwc(flag,b4nstflg))'*( temp(:,12,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsptbc(flag,b4nstflg))'*( temp(:,13,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsptwc(flag,b4nstflg))'*( temp(:,14,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsnw  (flag,b4nstflg))'*( temp(:,15,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) );
    grad_ngof_o(bwptflg,1) = -(ut.X2ptbc  (flag,bwpttflg))'*( temp(:, 3,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X2ptwc  (flag,bwpttflg))'*( temp(:, 4,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4sptbc (flag,bwpttflg))'*( temp(:, 8,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4sptwc (flag,bwpttflg))'*( temp(:, 9,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsptbc(flag,bwpttflg))'*( temp(:,13,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsptwc(flag,bwpttflg))'*( temp(:,14,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.Xngwptbc(flag,bwpttflg))'*( temp(:,16,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.Xngwptwc(flag,bwpttflg))'*( temp(:,17,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) );
    grad_ngof_o(bwftflg,1) = -(ut.X2ftbc  (flag,bwfttflg))'*( temp(:, 1,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X2ftwc  (flag,bwfttflg))'*( temp(:, 2,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4sftbc (flag,bwfttflg))'*( temp(:, 6,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4sftwc (flag,bwfttflg))'*( temp(:, 7,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsftbc(flag,bwfttflg))'*( temp(:,11,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsftwc(flag,bwfttflg))'*( temp(:,12,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.Xngwftbc(flag,bwfttflg))'*( temp(:,18,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.Xngwftwc(flag,bwfttflg))'*( temp(:,19,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) );
    grad_ngof_o(bwcflg,1)  = -(ut.X2ftwc  (flag,bwctflg ))'*( temp(:, 2,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X2ptwc  (flag,bwctflg ))'*( temp(:, 4,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4sftwc (flag,bwctflg ))'*( temp(:, 7,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4sptwc (flag,bwctflg ))'*( temp(:, 9,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsftwc(flag,bwctflg ))'*( temp(:,12,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsptwc(flag,bwctflg ))'*( temp(:,14,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.Xngwptwc(flag,bwctflg ))'*( temp(:,17,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.Xngwftwc(flag,bwctflg ))'*( temp(:,19,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) );

    grad_ngof_no=zeros(size(b));
    grad_ngof_no(1,1)       = -(ut.X2ftbc  (flag,cidx    ))'*( temp(:, 1,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X2ftwc  (flag,cidx    ))'*( temp(:, 2,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X2ptbc  (flag,cidx    ))'*( temp(:, 3,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X2ptwc  (flag,cidx    ))'*( temp(:, 4,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X2nw    (flag,cidx    ))'*( temp(:, 5,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4sftbc (flag,cidx    ))'*( temp(:, 6,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4sftwc (flag,cidx    ))'*( temp(:, 7,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4sptbc (flag,cidx    ))'*( temp(:, 8,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4sptwc (flag,cidx    ))'*( temp(:, 9,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4snw   (flag,cidx    ))'*( temp(:,10,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsftbc(flag,cidx    ))'*( temp(:,11,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsftwc(flag,cidx    ))'*( temp(:,12,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsptbc(flag,cidx    ))'*( temp(:,13,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsptwc(flag,cidx    ))'*( temp(:,14,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsnw  (flag,cidx    ))'*( temp(:,15,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.Xngwptbc(flag,cidx    ))'*( temp(:,16,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.Xngwptwc(flag,cidx    ))'*( temp(:,17,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.Xngwftbc(flag,cidx    ))'*( temp(:,18,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.Xngwftwc(flag,cidx    ))'*( temp(:,19,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no );
    grad_ngof_no(2,1)       = -(     (1-Beta).*gpdiff(flag, 1      ))'*( temp(:, 1,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (     (1-Beta).*gpdiff(flag, 2      ))'*( temp(:, 2,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (     (1-Beta).*gpdiff(flag, 3      ))'*( temp(:, 3,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (     (1-Beta).*gpdiff(flag, 4      ))'*( temp(:, 4,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (     (1-Beta).*gpdiff(flag, 5      ))'*( temp(:, 5,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (     (1-Beta).*gpdiff(flag, 6      ))'*( temp(:, 6,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (     (1-Beta).*gpdiff(flag, 7      ))'*( temp(:, 7,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (     (1-Beta).*gpdiff(flag, 8      ))'*( temp(:, 8,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (     (1-Beta).*gpdiff(flag, 9      ))'*( temp(:, 9,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (     (1-Beta).*gpdiff(flag,10      ))'*( temp(:,10,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (     (1-Beta).*gpdiff(flag,11      ))'*( temp(:,11,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (     (1-Beta).*gpdiff(flag,12      ))'*( temp(:,12,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (     (1-Beta).*gpdiff(flag,13      ))'*( temp(:,13,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (     (1-Beta).*gpdiff(flag,14      ))'*( temp(:,14,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (     (1-Beta).*gpdiff(flag,15      ))'*( temp(:,15,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (     (1-Beta).*gpdiff(flag,16      ))'*( temp(:,16,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (     (1-Beta).*gpdiff(flag,17      ))'*( temp(:,17,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (     (1-Beta).*gpdiff(flag,18      ))'*( temp(:,18,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (     (1-Beta).*gpdiff(flag,19      ))'*( temp(:,19,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no );
    grad_ngof_no(b2flg,1)   = -(ut.X2ftbc  (flag,b2tflg  ))'*( temp(:, 1,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X2ftwc  (flag,b2tflg  ))'*( temp(:, 2,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X2ptbc  (flag,b2tflg  ))'*( temp(:, 3,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X2ptwc  (flag,b2tflg  ))'*( temp(:, 4,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X2nw    (flag,b2tflg  ))'*( temp(:, 5,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no );
    grad_ngof_no(b4sflg,1)  = -(ut.X4sftbc (flag,b4stflg ))'*( temp(:, 6,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4sftwc (flag,b4stflg ))'*( temp(:, 7,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4sptbc (flag,b4stflg ))'*( temp(:, 8,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4sptwc (flag,b4stflg ))'*( temp(:, 9,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4snw   (flag,b4stflg ))'*( temp(:,10,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no );
    grad_ngof_no(b4nsflg,1) = -(ut.X4nsftbc(flag,b4nstflg))'*( temp(:,11,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsftwc(flag,b4nstflg))'*( temp(:,12,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsptbc(flag,b4nstflg))'*( temp(:,13,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsptwc(flag,b4nstflg))'*( temp(:,14,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsnw  (flag,b4nstflg))'*( temp(:,15,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no );
    grad_ngof_no(bwptflg,1) = -(ut.X2ptbc  (flag,bwpttflg))'*( temp(:, 3,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X2ptwc  (flag,bwpttflg))'*( temp(:, 4,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4sptbc (flag,bwpttflg))'*( temp(:, 8,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4sptwc (flag,bwpttflg))'*( temp(:, 9,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsptbc(flag,bwpttflg))'*( temp(:,13,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsptwc(flag,bwpttflg))'*( temp(:,14,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.Xngwptbc(flag,bwpttflg))'*( temp(:,16,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.Xngwptwc(flag,bwpttflg))'*( temp(:,17,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no );
    grad_ngof_no(bwftflg,1) = -(ut.X2ftbc  (flag,bwfttflg))'*( temp(:, 1,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X2ftwc  (flag,bwfttflg))'*( temp(:, 2,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4sftbc (flag,bwfttflg))'*( temp(:, 6,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4sftwc (flag,bwfttflg))'*( temp(:, 7,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsftbc(flag,bwfttflg))'*( temp(:,11,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsftwc(flag,bwfttflg))'*( temp(:,12,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.Xngwftbc(flag,bwfttflg))'*( temp(:,18,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.Xngwftwc(flag,bwfttflg))'*( temp(:,19,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no );
    grad_ngof_no(bwcflg,1)  = -(ut.X2ftwc  (flag,bwctflg ))'*( temp(:, 2,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X2ptwc  (flag,bwctflg ))'*( temp(:, 4,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4sftwc (flag,bwctflg ))'*( temp(:, 7,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4sptwc (flag,bwctflg ))'*( temp(:, 9,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsftwc(flag,bwctflg ))'*( temp(:,12,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsptwc(flag,bwctflg ))'*( temp(:,14,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.Xngwptwc(flag,bwctflg ))'*( temp(:,17,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.Xngwftwc(flag,bwctflg ))'*( temp(:,19,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no );

    grad_grno=zeros(size(b));
    grad_grno(1,1)       = -(ut.Xgwptbc (flag,cidx    ))'*( temp(:,16,3).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no ) - ...
                            (ut.Xgwftbc (flag,cidx    ))'*( temp(:,18,3).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no );
    grad_grno(2,1)       = 0;
    grad_grno(bwptflg,1) = -(ut.Xgwptbc (flag,bwpttflg))'*( temp(:,16,3).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no );
    grad_grno(bwftflg,1) = -(ut.Xgwftbc (flag,bwfttflg))'*( temp(:,18,3).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no );

    grad_grof_no=zeros(size(b));
    grad_grof_no(1,1)       = -(ut.Xgwptbc (flag,cidx    ))'*( temp(:,16,4).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.Xgwptwc (flag,cidx    ))'*( temp(:,17,4).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.Xgwftbc (flag,cidx    ))'*( temp(:,18,4).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.Xgwftwc (flag,cidx    ))'*( temp(:,19,4).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no );
    grad_grof_no(2,1)       = 0;
    grad_grof_no(bwptflg,1) = -(ut.Xgwptbc (flag,bwpttflg))'*( temp(:,16,4).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.Xgwptwc (flag,bwpttflg))'*( temp(:,17,4).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no );
    grad_grof_no(bwftflg,1) = -(ut.Xgwftbc (flag,bwfttflg))'*( temp(:,18,4).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.Xgwftwc (flag,bwfttflg))'*( temp(:,19,4).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no );
    grad_grof_no(bwcflg,1)  = -(ut.Xgwptwc (flag,bwctflg ))'*( temp(:,17,4).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.Xgwftwc (flag,bwctflg ))'*( temp(:,19,4).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no );

    grad_grof_o=zeros(size(b));
    grad_grof_o(1,1)       = -(ut.Xgwptbc (flag,cidx    ))'*( temp(:,16,4).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.Xgwptwc (flag,cidx    ))'*( temp(:,17,4).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.Xgwftbc (flag,cidx    ))'*( temp(:,18,4).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.Xgwftwc (flag,cidx    ))'*( temp(:,19,4).*(in_white_collar(flag)==1 | prev_WC(flag)==1) );
    grad_grof_o(2,1)       = 0;
    grad_grof_o(bwptflg,1) = -(ut.Xgwptbc (flag,bwpttflg))'*( temp(:,16,4).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.Xgwptwc (flag,bwpttflg))'*( temp(:,17,4).*(in_white_collar(flag)==1 | prev_WC(flag)==1) );
    grad_grof_o(bwftflg,1) = -(ut.Xgwftbc (flag,bwfttflg))'*( temp(:,18,4).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.Xgwftwc (flag,bwfttflg))'*( temp(:,19,4).*(in_white_collar(flag)==1 | prev_WC(flag)==1) );
    grad_grof_o(bwcflg,1)  = -(ut.Xgwptwc (flag,bwctflg ))'*( temp(:,17,4).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.Xgwftwc (flag,bwctflg ))'*( temp(:,19,4).*(in_white_collar(flag)==1 | prev_WC(flag)==1) );

    grad_no = grad_ngno+grad_grno+grad_ngof_no+grad_grof_no;
    grad_o  = grad_ngof_o+grad_grof_o;
    grad    = grad_no+grad_o;
end
